Criticality of a classical dimer model on the triangular lattice 
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We consider a classical interacting dimer model which interpolates between the square lattice case 
and the triangular lattice case by tuning a chemical potential in the diagonal bonds. The interaction 
energy simply corresponds to the number of plaquettes with parallel dimers. Using transfer matrix 
calculations, we find in the anisotropic triangular case a succession of different physical phases as the 
interaction strength is increased: a short range disordered liquid dimer phase at low interactions, 
then a critical phase similar to the one found for the square lattice, and finally a transition to an 
ordered columnar phase for large interactions. The existence of the critical phase is in contrast with 
the belief that criticality for dimer models is ascribed to bipartiteness. For the isotropic triangular 
case, we have indications that the system undergoes a first order phase transition to an ordered 
phase, without appearance of an intermediate critical phase. 
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I. INTRODUCTION 

Dimer models have regained a lot of interest in the con- 
densed matter physics community thanks to the pioneer- 
ing work of Rokhsar and Kivelson (RK) [l[ . In this orig- 
inal proposal, dimers represent singlets formed by pairs 
of spins 1/2 in Cuprate materials or frustrated antiferro- 
magnetic systems. This idea has been followed by many 
other proposals giving rise, at low energy, to such kind 
of effective models. To cite only a few connections to the 
physics of dimers, related again to magnetism, the lowest- 
energy configurations of fully frustrated Ising magnets 
can generally be mapped onto dimer configurations on 
the dual lattice 0, H| ■ Other specific (quantum) dimer 
models have recently been derived from a spin-orbital 
model describing LiNi02 [H , from the trimerized kagome 
antifcrromagnct Q, or for Heisenberg antiferromagnets 
under applied field in the pyrochlore lattice [f|. Lastly, 
we note that hard-core bosons or correlated fcrmions on 
frustrated lattices like the planar pyrochlore lattice can 
also be mapped onto dimer representations in the limit 
of large Coulomb repulsion 0, Q . 

Although the RK model is quantum mechanical, there 
is a special value of the parameters of the Hamiltonian 
for which the ground-state is an equal weight superposi- 
tion of all dimer tilings of the square lattice (RK point). 
Static properties at zero temperature can then be under- 
stood by the purely combinatorial problem of counting 
dimer tilings of the square lattice first solved in the early 
60's @. This system turns out to be a critical model, 
and more precisely a conformal field theory with central 
charge c = 1. Subsequently to the RK work, Moessner 
and Sondhi [Io| have shown that the same dimer model 
on the triangular lattice shows, instead of a single point 
with critical algebraic correlations, a disordered phase 
with short range dimer correlations - as expected from 
the known physics of its classical counterpart . From 
the interpretation originally given to the dimers, it is le- 
gitimate to call this phase a spin liquid. 

More recently, a study performed by Alet et al. [HJ 
on a square lattice classical dimer model showed the ex- 
istence of a Kostcrlitz-Thouless phase transition from a 



critical phase to a columnar ordered phase for the dimers. 
This transition is triggered by including an interaction 
term, consisting in the number of plaquettes doubly oc- 
cupied by parallel dimers. Subsequently, a quantum me- 
chanical dimer model can be built from this classical 
model proving that the single point showing criticality 
in the original RK Hamiltonian can be promoted to a 
whole critical phase (l3j . 

An important difference between the square and the 
triangular lattices for quantum dimer models lies in the 
degeneracy of the ground state at the RK point: it is 
finite for the triangular lattice with periodic boundary 
conditions (PBC) and scales exponentially with the lin- 
ear size for the square lattice with the same PBC (see 
Ref. [T3 and Ref. [l] respectively for details). This large 
degeneracy of the square lattice, related to its bipartite- 
ness, has been often designated as the responsible of the 
existence of critical correlations in this case. It is then 
natural to investigate the behavior of interacting dimer 
model on non-bipartite lattices interpolating between the 
square and the triangular lattices, and in particular try 
to elucidate the interplay between criticality and bipar- 
titeness. In the non-interacting case, it is found [Tl| that 
criticality disappears immediately with the introduction 
of non-bipartite dimers. 

In this paper we study an extension of the classical 
dimer model with nearest-neighbor interactions already 
studied elsewhere [lU, [l3| to the case of an anisotropic tri- 
angular lattice. Although a limiting case of this model is 
the square lattice, in general the lattice is non-bipartite. 
In Sec. II we start by reviewing the results on the square 
lattice, and recalling the principles of an associated field 
theory used to describe the critical phase. In Sec. Ill, 
we investigate the properties of the anisotropic model 
that interpolates between the square and triangular lat- 
tices by means of transfer matrix calculations. There are 
two independent paths for the interpolation: the intro- 
duction of a fugacity for diagonal bonds but also of the 
introduction of interactions between parallel dimers on 
diamond-like plaquettes (see Fig. [1}. We find for a cer- 
tain range of parameters the existence of a critical phase. 
We determine its location in the phase diagram as well 
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as transitions from it to either liquid or ordered phases. 
In Sec. IV, we find that the isotropic model on the tri- 
angular lattice, does not possess a critical phase for any 
interaction strength but could display an ordered phase 
for sufficiently large interactions. Finally, Sec. V con- 
tains a discussion of the various results and conclusions. 




(a) (b) 



FIG. 1: (Color online) Dimer models a) on the square lat- 
tice, with dimer interactions u on square plaquettes; b) on an 
anisotropic triangular lattice, with dimer interactions u and 
v on square and diamond-like plaquettes respectively and a 
fugacity zd for dimers on diagonal links. 



II. CRITIC ALITY OF THE DIMER MODEL ON 
THE BIPARTITE SQUARE LATTICE 



Let us start by reviewing the main results of the square 
lattice dimer model. In this paper, we will only consider 
fully-packed dimers with a hard-core constraint: there 
must be one and only one dimer coming out of each 
site. The non-interacting model can be solved exactly 
on the square lattice by expressing the partition function 
of the system as a Pfaffian, which enables to obtain an- 
alytic expressions for the free energ y a nd correlators, on 
either finite or infinite systems [3, [l4[. On the square 
lattice, dimers are found to be in a critical phase: dimer- 
dimer correlations decay algebraically as 1/r 2 with rel- 
ative distance r, whereas the monomer-monomer (defin- 
ing a monomer as a site with no dimer) correlations be- 
have as 1/ y/r. Classical dimers on the square lattice have 
also been studied within a finite-temperature model with 
nearest neighbors interactions, i.e. an interaction —u 
between parallel dimers on the same plaquette (u > 
corresponding to attractive interactions) [12|]- Here the 
model is not intcgrable anymore and it was found nu- 
merically that the system is critical down to a finite 
temperature u co \/T = 1.54(3), where it gives rise to a 
low-tempcraturc columnar phase with dimers aligned in 
columns. Monomer-monomer correlations decay to at 
long distances, either algebraically (in the critical phase) 
or exponentially (in the columnar phase). From now on, 
unless it is specified otherwise, we set T = 1 and let 
the different dimer-dimer interactions (u for the purely 
square lattice) vary. Nevertheless, we still refer to high-T 
(respectively low-T) regions characterized by small (resp. 



large) values of the involved couplings. 

The criticality encountered in the high-T phase can be 
understood in the framework of a height field theory. For 
each dimer configuration on the square lattice, a scalar 
height field with quantized values (multiples of j) can 
be defined microscopically on the dual lattice [l5| f. The 
spatial variations of the height field between neighboring 
sites of the dual lattice are entirely determined by the 
presence of dimers between them, thanks to the lattice 
bipartiteness. The long-wavelength modes of this height 
correspond to a coarse-grained height field $(r) defined 
in continuum space, and the physics of the model is cap- 
tured by the action [l2|, EH : 

S[$\= y rf 2 r ff 7r|V$(r)| 2 + Fcos(27rp$(r)) (2.1) 

The cosine term of this action is a locking potential that 
favors p flat configurations (p = 4 for the square lattice), 
corresponding microscopically to the columnar configu- 
rations. The | V*!? | 2 term, corresponding to the cost of 
fluctuations around these flat configurations, account for 
the entropy of dimer coverings. In the non-interacting 
case, the value of the stiffness is fixed by the exact re- 
sults of Ref. @, [3 to be g = \ (see also Ref. [13). The 
renormalization of the stiffness constant accounts for the 
role of interaction: g increases with u. 

As long as the cosine operator in Eq. (|2.ip is irrele- 
vant in the renormalization group sense, this action de- 
fines a conformal field theory with central charge c = 1 
and the system is critical. The dimer operator (which 
gives the local dimer density) is composed of a cos(27r$) 
and gradient terms (l8j . and has a scaling dimension 
(ii^ = 2^ , corresponding to half of the exponent of 
dimer-dimer correlations. Similarly, the locking poten- 
tial has a dimension d„ n = = -, and becomes rel- 
evant when g > 4. Within the standard Coulomb gas 
description [19( of this theory, we can define operators 
corresponding to the insertion of a particle of electric 
and magnetic charges e and m respectively. In general, 
the dimension of such an electromagnetic vertex opera- 
tor is d e . m = + In the height model discussed 
here, the operator with (e,m) = (1,0) is the dimer op- 
erator already mentioned; similarly, the monomer oper- 
ator corresponds to inserting a magnetic particle with 
(e,ra) = (0, ±1) (the sign of the charge depends on the 
sublattice where the monomer is inserted). In the criti- 
cal phase, dimer-dimer and monomer-monomer correla- 
tors decay as power laws of the distance and the decay 
exponents are 2(^0 and 2do.i respectively. Going back to 
the height picture, the insertion of magnetic charges (for 
instance monomers or links breaking the bipartiteness) 
corresponds to dislocations and, as we show below, can 
be treated as perturbations in the coarse-grained field 
theory. 

Prior to the presentation of our strategy to attack the 



3 



dimcr problem, wc first introduce the numerical tech- 
niques used here and illustrate them on the example 
of the square lattice interacting dimer model. In this 
study we use a Transfer Matrix (TM) approach to de- 
termine the domains of existence of the different phases 
of the dimer model, computed on a torus of longitudi- 
nal and transverse sizes L T and L respectively. The TM 
techniques allow to compute exponents of critical dimer- 
dimer and monomer- monomer correlations and also cor- 
relations themselves (see Appendix) , on systems with L T 
large enough to be considered infinite. 

For a critical phase such as the one encountered in 
the high-T region of the square lattice dimer model, we 
can compute the decay exponents of correlation func- 
tions in two ways: either by a direct inspection of the 
real-space decay of correlation functions (see the descrip- 
tion in Sec|3]of Appendix), or by the leading eigenvalues 
of the TM (see SecQ]of Appendix). Both methods are 
guided by a conformal field theory (CFT) analysis of the 
critical phase, and the latter is more precise to deter- 
mine exponents as it allows to use efficiently topologi- 
cal and translation symmetries of the TM (see SecQ] in 
Appendix). The calculations of correlation functions in 
real space offer however the advantage of providing more 
physical informations (for example about the correlation 
length) when the system is not critical. Another impor- 
tant quantity, namely the central charge c of the CFT, 
can also be extracted from the finite-size scaling of the 
TM largest eigenvalue when the system is critical (see 
SccQ]of Appendix). 

We illustrate the validity of this numerical approach in 
Fig. [2] where the critical phase is evidenced by the c = 1 
plateau of the central charge (estimated using TM largest 
eigenvalues for system widths up to L = 14). The scal- 
ing dimensions of dimer d\fi and monomer c?o,i operators 
(obtained from sublcading eigenvalues) are also displayed 
according to the discussion above, the transition to the 
columnar phase is characterized by a value g = 4 of the 
stiffness, thus a value di.o = 1/8. This criterion is used 
to estimate the transition temperature (estimates of do,i 
can also be used in principle, but they are more sensi- 
tive to finite-size effects in that temperature range). The 
results of Fig. [2] lead to an estimate u nr >i = 1.563(7), in 
good agreement with previous results [lj, [H| . 



III. CRITICAL PHASE IN A DIMER MODEL 
INTERPOLATING THE SQUARE AND 
TRIANGULAR LATTICES 



In this section, we investigate whether the bipartite- 
ness of the square lattice is a condition for the existence 
of the critical phase, by building a model defined on a 
lattice interpolating continuously between the square and 
the triangular lattices. 
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FIG. 2: (Color online) Estimates and extrapolations from dif- 
ferent estimates (see Appendix for details and notations) of 
the central charge c and of the scaling dimensions di,o and 
do,i in the dimer model on the square lattice as a function of 
interaction strength it. 



Definition of a model with lattice and interaction 
anisotropics - For commodity, we represent the tri- 
angular lattice as a deformed square lattice with bonds 
in the x, r and t — x direction (sec Fig. [T]and Appendix). 
The interpolation between both lattices is made by as- 
signing a fugacity zd to the diagonal bonds. Since the 
elementary plaquettes of the triangular lattice are more 
numerous than those of the square, we need to define 
another parameter v characterizing the interaction on 
diamond-like plaquettes (i.e. non-square four-site pla- 
quettes). Wc finally keep the notation u for the inter- 
action between parallel dimers on a square plaquettc. 
These elementary energy scales in the problem are il- 
lustrated in Fig. [TJ For zd = and v = 0, we have the 
classical dimer model on the square lattice discussed in 
Sec. HH while the pure triangular dimer model (that will 
be discussed in Sec. IIV|) corresponds to both zd = 1 and 
v/u = 1. To describe the system in function of these two 
anisotropy parameters, we take the approach of varying 
cither zd or v/u while keeping the other parameter con- 
stant. The different paths in the (zd,v/u) space where 
phase diagrams were computed numerically are repre- 
sented in Fig. [3J 

Description of expected phases and transitions — In 
the non-interacting limit (u = v = 0) of the present 
model, the system is disordered with exponentially de- 
caying correlations for any fugacity Zjj > 0: criticality 
is destroyed by an infinitesimal proportion of diagonal 
bonds [llj . This can be accounted for in the height field 
theory by adding a relevant term in the action Eq. (|2.ip . 
which can be understood by noticing that adding a sin- 
gle diagonal dimer on the square lattice corresponds to 
inserting two monomers on the same sublattice, i.e. two 
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FIG. 3: (Color online) Different parameter-space paths (A to 
E) taken in the TM calculations in section Hill 



magnetic charges of the same sign [20[ . The correspond- 
ing m = 2 vertex operator is written in terms of the 
field 8, dual to the height field <J>, and leads to the new 
effective action: 

S[$} = J d 2 r7r, 9 |V$(r)| 2 + ^cos(87r$(r))-Acos(2e(r)) 

(3.1) 

In the absence of interactions, this new perturbing field 
(the cos(20) term) corresponds to the mass term of the 
free Majorana doublet arising in the large scale re gim e 
within the Pfaffian description of the dimer model [ll| . 
For a generic value of the stiffness, the scaling dimensions 
of the two perturbing terms are e?4.o = | and do, 2 = 2g 
respectively. The fact that the introduction of diagonal 
bonds is indeed relevant in the non-interacting limit is 
seen from the value of the stiffness in this case g = h 
(since do, 2 = 1 < 2). However, the key point of our 
study is that the cos(20) term becomes irrelevant for 
g > 1 whereas the locking potential is relevant only for 
g > 4. As we can tune the stiffness constant by modifying 
the interaction strength between parallel dimers (u and 
v terms), we are hoping to reach the 1 < g < 4 window, 
i.e. a critical phase on the triangular lattice, within the 
interacting dimer model defined above. The predicted 
phase diagram is therefore as follows: first, a high-T liq- 
uid phase (for 1/2 < g < 1 in the unperturbed model), 
then a transition to a critical phase (1 < g < 4) at inter- 
mediate couplings and finally a low-T columnar ordered 
phase (for g > 4). The dimer d^o and monomer do,i scal- 
ing exponents are predicted to be both equal to 1/2 at 
the transition from the liquid to the critical phase, and 
1/8 and 2 respectively at the critical-columnar transition. 
In both cases, transitions out of the critical phase are ex- 
pected to be of Kosterlitz-Thouless type. This analysis 
holds only in the perturbative regime for the diagonal 
dimer fugacity Zd *C 1, and we now address the question 
whether this scenario is realized for an arbitrary lattice 
anisotropy < Zr> < 1 by means of numerical TM calcu- 
lations. 



Perturbation of the square lattice model by diagonal 
bonds — We first consider the model with no additional 
interactions v = and turn on the diagonal dimer fugac- 
ity (path A in Fig. [3]). Our numerical estimates of the 
central charge c and exponents d\$ and do.i are displayed 
in Fig. [5] as a function of the coupling strength u for the 
specific value zr> = 0.4. We clearly observe the emer- 
gence of a c = 1 plateau witnessing a critical phase, for 
a wide range of couplings. The points u\ and u 2 where 
di,o = 1/2 and do,i = 1/2 correspond roughly to the 
high-T limit of the c = 1 plateau, as predicted above; 
from these criteria, the transition between the critical 
and liquid phase is located at u* = (ui+u 2 )/2 = 0.80(5), 
the error bar being estimated by (u 2 — u{)/2. Similarly, 
the low-T end of the plateau at u co i = 1.566(8), corre- 
sponding to the entrance into the columnar phase, coin- 
cides with the criterion d^o = 1/8, confirming the analy- 
sis above. We use estimates of difi to locate the columnar 
transition, since as on the square lattice estimates of the 
monomer exponent <io,i are more affected by finite-size 
effects. 
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FIG. 4: (Color online) Central charge c, estimates and extrap- 
olations of the dimer exponent di t o and monomer exponent 
do,i, as a function of it for a diagonal dimer fugacity zo — 0.4 
and v = 0. 



Repeating the same analysis for different values of the 
diagonal dimer fugacity, we obtain, as a function of zr> 
and u, the phase diagram of Fig. [5l The main features 
are: (i) the transition temperature to the columnar phase 
is essentially not affected by the diagonal dimer fugacity 
and remains close to the value obtained in the square 
lattice; (ii) as found in Ref . [Ill, the system with no inter- 
action u = is a gapped liquid irrespective of the diag- 
onal dimer fugacity; (Hi) the dimer model (with v = 0) 
comprises a critical phase even on the isotropic triangu- 
lar lattice (zr> = 1), where it extends from u* = 1.08(5) 
up to u co i = 1.575(10). 

These results are in full agreement with the predictions 
of the field-theoretical analysis developed above. We em- 
phasize also here that the lattice bipartitcness is not a 
necessary condition to have criticality in a dimer model, 
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and that this condition should rather be replaced by a 
condition on the existence of a stiffness constant window 
1 < g < 4 in an unperturbed model. 
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FIG. 5: (Color online) Phase diagram for v — (line A in 
Fig. I3J) as a function of diagonal dimer fugacity zd and squarc- 
plaquettes interaction u. Note that on the zd = line, the 
small u region also corresponds to a critical phase (dashed 
line) . 



Introduction of a finite interaction anisotropy v/u - 
We now introduce the interaction term v on the diamond- 
like plaquettes. The u and v terms are in competition, 
since they tend to favor different dimer orderings. When 
v < u, the square columnar configurations will eventu- 
ally dominate at low enough T, but the presence of the 
frustrating v interactions will shift down the transition 
temperature to the ordered phase. As the algebraic cor- 
relations of the critical phase also correspond to a square- 
like columnar ordering at lower temperatures, we also 
expect the extent of the critical phase to shrink down as 
the strength of the v interactions is increased. From the 
point of view of the effective action Eq. (|3.1[) . we expect 
g to decrease with the v frustrating interactions. 

For a sufficiently small v/u parameter, the behav- 
ior of the system is found to be identical to the case 
v = 0, with a critical phase, evidenced by a c = 1 
plateau between the dimer liquid and the ordered phase. 
This is illustrated in Fig. [6] for anisotropy parameters 
{zd, v/u)=(0A, 0.2). As expected, the transition towards 
the ordered phase is slightly shifted towards lower tem- 
peratures (ticoi = 1.624(8) for v/u = 0.2 while u co i = 
1.566(8) for v = 0). The high-T boundary of the critical 
phase is more strongly affected by the presence of v in- 
teractions: for the same lattice anisotropy, u* = 0.92(6) 
for v/u = 0.2 to be compared with 0.80(5) for v = 0. 

Using the same procedure as before, the transition 
points u* and u co \ were determined in the regions of the 
phase diagram corresponding to paths B,C,D (with re- 
spectively v/u — 0.2, 0.4 and 0.5), and E (for zd — 1) and 
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FIG. 6: (Color online) Central charge and exponents (es- 
timates and extrapolations) for an interaction anisotropy 
v/u = 0.2 and diagonal dimer fugacity zd = 0.4. 
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FIG. 7: (Color online) Phase boundaries u* and func- 
tion of either zd [(a), (b) and (c)] or v/u [(d)], corresponding 
respectively to the paths B, C, D and E in Fig. [3] Error bars 
on u co i are smaller than symbol size. 



the corresponding phase diagrams are reported in Fig. [7] 
As anticipated, we clearly see that the effect of adding 
the interaction v is essentially to shift both boundaries of 
the critical phase towards lower temperatures: in other 
words, the frustration due to the v terms has a net ten- 
dency to stabilize the dimer liquid and to destabilize the 
ordered (columnar) or quasi-ordered (critical) phase. For 
a given fugacity Zd, the extent of the critical phase indeed 
decreases as v/u increases. The transition between the 
dimer liquid and the critical phase turns out to be more 
affected by v/u than the transition towards the columnar 
phase. We can interpret this with the following rough ar- 
gument: v interactions effectively increase the weight of 
diagonal bonds in the dominant configurations, enhanc- 



6 



ing the perturbation caused by the diagonal bonds and 
therefore decreasing the extent of the critical phase from 
the high-T direction. 

For large values of v/u, a fair determination of phase 
boundaries becomes more difficult as finite-size effects are 
getting more important. For instance, the determination 
of u* with the dimer and monomer exponents criteria 
give values increasingly far from each other, resulting in 
larger error bars in the phase diagrams of Fig. [7t and 
d. In practice for v/u > 0.5, both our numerical results 
and the way of analyzing them break down due to uncon- 
trolled finite size effects. This is actually not a surprise, 
as the CFT-guided analysis of the TM results (see Sec. 
[5] in Appendix) relies on the existence of an unperturbed 
critical window in the model, which is no longer present 
for large v interactions. For v/u > 0.5, we also find that 
the second largest TM eigenvalue is no more found in the 
q = it symmetry sector in the low-T phase. This indicates 
a change in nature of the lowest-energy excitations and 
that the scheme along which criticality was previously 
understood is not valid anymore. From these arguments 
and even if our finite-size results become less reliable, it 
is clear that they are no longer compatible with the ex- 
istence of an extended critical phase for v/u > 0.5: this 
is exemplified by the absence of a plateau in the numer- 
ical estimate of the central charge. Finally, we note that 
in parallel to the shrinking of the critical phase observed 
when increasing v/u, the columnar transition seems to be 
progressively shifted towards the T — limit as the inter- 
action anisotropy parameter approaches 1 (see Fig. Hi). 
This is again in agreement with our previous qualitative 
arguments on the effect of v interactions. 



IV. DIMER ORDERING ON THE ISOTROPIC 
TRIANGULAR LATTICE 



We now turn to the more complex situation of the 
isotropic triangular model, where Zd = 1 and v = u 
(in this section, the isotropic interaction will only be de- 
noted by u). The high-T limit is in this case well under- 
stood: the system is in a dimer liquid phase, with expo- 
nentially decaying dimcr-dimer and monomer- monomer 
correlation functions. This has been shown exactly for 
u = v = in Ref. [ill , and we could confirm it for 
small values of the interaction strength. For intermediate 
u = v values, as expected from the results of Sec. IIIII for 
v/u > 0.5, we find no evidence for a critical phase in our 
TM calculations (in particular, noc=l plateau). 

Whereas we always found a low-T columnar ordering 
up to now, the low-T behavior of the isotropic trian- 
gular model reveals more complex. Indeed, from the 
isotropy of the interactions and since all dimer fugaci- 
ties are equal, it appears that configurations that min- 
imize the energy are more numerous. At first glance, 



12(= 6x2) columnar ordered ground-state configura- 
tions emerge (6 for the number of ways to put a dimer 
on a link connected to a given site and 2 for the possi- 
ble directions of a dimer column, given an orientation of 
dimers). However, as first remarked in Ref. [ToL the num- 
ber of configurations that minimize the energy is much 
larger as one can easily create zero-energy defects by just 
translating a line of dimers or flipping all dimers along a 
column. These two types of moves from a given columnar 
ground-state are represented in Fig.[H] The almost exten- 
sive ground-state degeneracy generated by these costless 
defects clearly changes the picture for the low-T behav- 
ior of the isotropic triangular dimer model. Indeed, at 
T = 0, strictly speaking, the system is no longer ordered 
with respect to a local order parameter such as the colum- 
nar one defined for the square lattice. However, it is still 
possible that at finite T, the thermal fluctuations select 
a specific ordering pattern: this would be an illustration 
of the "order by disorder" effect [MUl]. 



(a) (b) 




FIG. 8: (Color online) Different types of zero-energy defects 
on the (isotropic) triangular lattice dimer model, from a spe- 
cific columnar ordered ground-state: (a) line-shifting modes; 
(b) column-flipping modes. 



In fact, for the quantum dimer model on the triangular 
lattice, Moessner and Sondhi found in perturbation the- 
ory (with respect to the quantum kinetic terms) that in- 
deed (quantum) fluctuations select the 12 ordered colum- 
nar states 0. Noticing that in our classical model, ther- 
mal fluctuations play the same role as quantum kinetic 
fluctuations in the quantum dimer model [l2l ] since both 
essentially count the number of flippablc plaquettes in 
a given state, it is likely that thermal fluctuations (in- 
stead of quantum fluctuations) trigger a similar order- 
by-disorder scenario here. From general arguments [22| . 
we expect that if it is indeed the case, the transition 
temperature to the columnar phase should be quite low. 
It could also happen that, for our specific model, this 
transition never occurs and the systems stays in a liquid 
phase down to T = 0. 
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FIG. 9: (Color online) Connected dimer-dimer correlations as 
a function of the distance r T along the main axis of a torus 
of dimensions (L T ,L) = (60,6), for interaction strengths it 
between 1 and 5. Correlations of vertical (respectively hori- 
zontal) dimers are represented by up (resp. down) triangles. 
For clarity, we only displayed horizontal (resp. vertical) dimer 
correlations at odd (resp. even) distances. Inset: estimated 
correlation length £ as a function of u. 



peak but already note at this stage that the clear nar- 
rowing of the specific heat peak with increasing system 
size, conjugated with the absence of a power-law envelop, 
is very suggestive of a first- order phase transition. 




FIG. 10: (Color online) Specific heat per site c v = L X L d ^p 
as a function of u for systems of different widths (see Ap- 
pendix for details of computation). 



We now try to settle this issue with the help of numer- 
ical TM calculations. As the critical phase disappeared, 
we can no longer use the CFT-based analysis of the previ- 
ous sections and we have to resort to more standard ther- 
modynamical means of detecting the hypothetic transi- 
tion to a low-T ordered phase. We first computed dimer- 
dimer correlations by TM iterations (see Appendix). In 
Fig. connected correlation functions are displayed as 
a function of dimer-dimer distance r T for the example of 
a (L T , L) = (60, 6) torus for interaction strength u up to 
5. As exemplified by the log-linear scale, they are clearly 
short-ranged with an exponential decay (plus an oscillat- 
ing part depending on the parity of r T ). We obtain an 
estimate of the correlation length £ by fitting the correla- 
tions to an exponential decay exp(— r T /£). The resulting 
curve £(w), displayed in the inset of Fig. O shows a well- 
pronounced peak at a finite value of u, which could corre- 
spond to an ordering temperature. We also find that the 
peak positioned at u = 1.9(1) for the system size L = 6, 
shifts towards higher interaction strength u = 2.22(5) for 
L = 8. 

Another thcrmodynamical insight is given by the spe- 
cific heat per site c v = L L d{E)/dT , which is also ac- 
cessible to the TM calculations (see Sec. [2] of Appendix 
for details). Fig. [TU] displays the specific heat as a func- 
tion of u for different system widths L. For each value 
of L, the specific heat peaks at a temperature close to 
the ones of the £(u) curves. As L increases, the peak 
in c v (u) sharpens and shifts to lower temperatures. We 
discuss below the finite-size scaling of the position of the 



Having possible signs of a phase transition, we now 
try to find an order parameter for the low-T phase. 
As already discussed, this a not a simple task because 
of the zero-energy modes that are responsible for the 
large ground-state degeneracy. The specific geometry 
of the TM is of help here to characterize the appear- 
ance of long-range order. Consider a very long cylin- 
der L T ^> L ^> 1 with PBC in the small L direction. 
With this geometry, the ground-states having dimers per- 
pendicular to the long direction (horizontal dimers) are 
much more numerous than ground-states with only ver- 
tical or diagonal dimers. Indeed, there are 2 Lt horizontal 
ground-states (corresponding to the line-shifting modes, 
see Fig. |8^). The same line-shifting modes provide only 
a much smaller (2 L ) degeneracy for vertical or diago- 
nal ground-states. The other family of low-energy modes 
(column- flipping modes) give a small 0{2 L / 2 ) de gener- 
acy for the horizontal ground-states (see Fig. [8)3), and 

0{2 L -/' 2 ) dc generacy for both vertical and diagonal 
ground-states. Consequently, in the limit of large L T , 
we expect horizontal line-shifting modes to predominate. 

If there is long-range order in the system, the previ- 
ous analysis indicates that the long cylinder geometry 
induces a preferential ordering in the horizontal dimers 
(for entropic reasons). Consequently, we consider the av- 
erage occupation of horizontal bonds on a given line of 
the lattice as a good indicator of a possible phase tran- 
sition. At infinite temperature, occupations of all bonds 
are all likely and equal to 1/6. In contrast, from the ar- 
guments above, at very low T, the probability for a hori- 
zontal bond to be occupied on an infinitely long cylinder 
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is P_ = 1/2. Defining 



P_ - 1/6 



1/2 



= 3P_ - -, 
1/6 2' 



we expect the "order parameter" (m) to vanish in the 
low-coupling limit u = and to saturate to 1 at large 
enough coupling. Even if. strictly speaking, this above 
argument is not rigorous (because L T and L are both 
finite in our computations), we expect (m) to reflect a 
true physical behavior, (m) is shown in Fig.[TT]for lattice 
widths L = 4, 6, 8 and we clearly observe the predicted 
behavior (saturation to at small u and to 1 at large 
couplings). As for the specific heat data, we observe that 
the shape of (m) is strongly influenced by the system 
size, with similar trends: as L increases, the temperature 
region of ordering (where d(m) / du is maximal) is getting 
narrower and is shifted towards lower temperatures. 



<m> 



% — & — ©- 




3.5 



FIG. 11: (Color online) Proposed order parameter (m) for the 
long cylinder geometry (see text) versus coupling strength u 
for different lattice sizes. Inset: finite-size transition temper- 
atures and T)? (maxima of c v and d{m)/du respectively) 
versus 1/L. 



In order to have a better understanding of this ordering 
phenomenon, we have computed the order parameter (m) 
away from the purely isotropic case, for a diagonal dimer 
fugacity zjj = 1 but for values of v/u lower than 1. The 
evolution of (m) with u is shown on Fig |12l for different 
interaction anisotropies v/u = 0.4; 0.6; 0.8; 1. As soon as 
u is larger than v, the mean occupation of horizontal links 
P_ is 1/4 at zero temperature, which is consistent with 
the fact that, as in the purely square columnar state, only 
4 configurations have the minimal energy. Nevertheless, 
the ordering occurring when u increases is characterized 
by larger values of (m) in an intermediate temperature 
range. As v/u increases, the maximal value of P_ gets 
closer to 1/2, which is the value found in the ordered 
phase discussed previously in the isotropic case. Notice 
that the extent of the temperature range where P_ ap- 
proaches 1/2 increases with the length L T of the system 



(see v/u = 0.6 curves in the Fig. I12p. which supports the 
scenario of order by disorder: the P_ = 1/2 asymptotic 
value can be understood by the number 0(2 Lt ) of lowest- 
energy defects (which are lowest-energy configurations in 
the isotropic case). With interaction anisotropy, these 
defects are allowed only at finite temperatures, and pro- 
liferate down to lower temperatures as u — v decreases. 
This clearly shows that the ordered phase in the isotropic 
triangular lattice model differs both qualitatively and 
quantitatively (by the number of lowest-energy config- 
urations) from that of its anisotropic version, and is a 
good example of the "order by disorder" scenario. 




v/u=0.4 
v/u=0.6 

v/u=0.6 (L T =150) 

v/u=0.6 (L=200) 

v/u=0.8 
v=u 



3 



U 



FIG. 12: (Color online) Average occupation P_ of horizontal 
bonds computed on a row of an isotropic triangular lattice of 
dimensions L = 6 and L T — 100 (if not specified) , for different 
interaction anisotropy parameters — = 0.4; 0.6; 0.8; 1, as a 
function of interaction strength. 



At this stage, we have three indications (from the 
correlation length, the specific heat and (m)) that the 
isotropic triangular interacting dimer model could order 
at low temperatures. In all cases, the transition tempera- 
ture was seen to decrease quite consequently with system 
size. To check whether the system could order at finite 
T in the thermodynamic limit, we perform a finite-size 
scaling of these effective transition temperatures. In the 
insert of Fig. 111! the finite-size temperature transitions 
T*(L) and T%(L) corresponding to respectively the max- 
ima of the specific heat and of d(m)/du, are plotted as 
a function of inverse transverse system size 1/L. A pre- 
cise finite-size scaling form is difficult to determine (due 
noticeably to large error bars), but all reasonable finite 
size dependence (e.g. linear or quadratic in 1/L) lead to 
a finite value of T c in the thermodynamic limit. A rough 
estimate can be made with the help of linear interpola- 
tion (dashed line on the figure) and gives T c = 0.2 ± 0.05 
(in units of u). 

Our numerical results therefore seem to be consistent 
with a finite temperature ordering of the isotropic lat- 
tice model, probably triggered by an order by disorder 
mechanism |2ll . l22j ]. The sharp behavior of both spe- 
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cific heat and "order parameter" at the transition (see 
Fig. [TTj) , as well as the absence of any criticality behav- 
ior in the central charge, suggests that this transition is 
first order. We finish by noting that the samples used 
in the computations are of relatively moderate size, and 
that it is still possible that the extrapolated transition 
temperature actually vanishes with larger samples avail- 
able. To have a better understanding of this transition 
with the presently available system sizes, one could also 
study this model on samples with other geometries and 
adapted order parameters, and check if extrapolations 
of finite-size temperatures on different geometries give 
the same T c in the thermodynamic limit. Another in- 
sight could be given by Monte Carlo simulations for this 
model. We expect however these simulations to be dif- 
ficult as the presence of the low-energy modes (and the 
corresponding large degeneracy of the ground-state) will 
certainly induce ergodicity and freezing problems in the 
Monte Carlo process [23| . Such investigations are beyond 
the scope of the present paper. 



V. CONCLUSIONS 



To summarize, we have constructed a simple classi- 
cal interacting dimer model on a lattice that interpolates 
between square and triangular lattices. This is of par- 
ticular interest since the two limiting models singularly 
exhibit very different behaviors at infinite temperature 
(critical and short-ranged phases for the square and tri- 
angular lattices respectively). Since the topology of the 
triangular lattice can be simply obtained from the square 
lattice by adding one extra diagonal bond on each square 
plaquette, we have introduced a fugacity parameter for 
these extra bonds. Similarly, the extra bonds leading 
to new local interactions in diamond-like plaqucttes, we 
have also considered dimer interactions on these plaquc- 
ttes, differing from those on the square plaquettes. 

This anisotropic dimer model has been investigated in 
great details using simple considerations as well as nu- 
merical transfer matrix techniques (on strips with up to 
14 sites wide) supplemented by predictions from confor- 
mal field theory. A very rich phase diagram has been 
obtained with, in particular, a novel intermediate behav- 
ior where the critical phase is now restricted to a finite 
intermediate temperature range and does not extend, as 
for the pure square lattice model, up to infinite temper- 
ature. At high-temperature, the critical phase is hence 
replaced by a liquid dimer phase as in the isotropic tri- 
angular lattice. It is of interest that such a behavior 
appears immediately for arbitrary small fugacity of di- 
agonal bonds so that one can view the typical behavior 
of the square lattice as a "singular" limit. However, we 
note that the critical phase itself survives in an extended 
vicinity of the square lattice model (although its exten- 




v/u 

1 

FIG. 13: Schematic projective phase diagram of the dimer 
model on the triangular lattice (zd = 1, line E of Fig. |3J, in 
terms of the interaction u and interaction anisotropy v/u. 



sion in temperature jumps abruptly) and this stability 
is therefore not related to the bipartiteness of the lat- 
tice. While the critical phase disappears progressively 
with increasing interaction anisotropy, the ordered phase 
survives for any range of anisotropy, as well as the high-T 
liquid dimer. Criticality in this model hence requires to 
have sufficiently anisotropic interactions. These findings 
arc reproduced in the schematic phase diagram of Fig. 1131 
for the case of isotropic dimer fugacity zjj = 1. Lastly, 
we have devoted special attention to the limiting case of 
the isotropic dimer model on the triangular lattice where 
we found evidence of a direct (likely first-order) transi- 
tion between the high-T liquid phase and a low-T ordered 
phase triggered by an order-by-disorder mechanism. 
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TRANSFER MATRIX ANALYSIS OF THE 
ANISOTROPIC DIMER MODEL 



Other works [TJ, [13| already discussed the construction 
of the TM for an interacting dimer model. We give a few 
technical details related to the triangular lattice and to 
the observables studied in the body of the text. 
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L 


dim(T (q = 0)) 


dim(T (q = it)) 


l ■ / mn 1 r\\\ 

dim(T (q = 0)) 


8 


424 


421 


410 


10 


2980 


2929 


2954 


12 


22218 


22207 


22150 


14 


170980 


170665 


170822 



TABLE I: Dimensions of the TM restricted by translation 
invariance to symmetry sectors of wave vectors q = and 
q = 7r, and by topological invariance to even and odd sectors, 
for system widths L up to 14. 



1. ID and 2D dimer configurations 



The TM allows to treat a 2D Lx L T system with peri- 
odicity in the transverse direction (x) as a succession of 
L T ID rows of lengths L. In a row at abscissa r (the r axis 
being the direction of propagation of the TM), a dimer 
configuration C is specified by the dimer occupation of 
horizontal (oriented along e^) links at abscissa r and of 
vertical and diagonal links connecting sites at abscissa r 
to sites at r+1, at relative positions e T and e T e x respec- 
tively. The number of such configurations, which is the 
size of the TM, is 3 L (This numbe r can be evaluated as 
Tr(M L ) where M is a one-dimensionnal transfer matrix 
defined between the 4 dimer configurations allowed on a 
triangle, taking into account the hard-core constraint). 

The topology of the triangular lattice is such that on 
a i-wide cylinder, the parity of the number of vertical 
or diagonal dimers linking two successive rows is con- 
served along the r axis (assuming that L is even). This 
defines two topological sectors [even and odd), and the 
TM is block-diagonal with corresponding blocks T e and 
T° respectively. We also use the invariance by transla- 
tion along the x axis to reduce the size of T. Translation 
invariance allows to consider, instead of 3 L ID configu- 
rations, only their representatives in a symmetry sector 
(defined by the transverse wave vector q), from which 
other configurations are obtained by translations along 
e x . Note that the TM elements take into account not 
only representatives but also all configurations connected 
to a given representative. These two invariances allow 
to reduce considerably the TM size in various symme- 
try and topological sectors. We give in Table[T]the size of 
the different sectors useful in the following for L up to 14. 
Values of non-zero TM elements, i.e. between compati- 
ble ID-configurations, depend on the number of doubly 
occupied plaquettes of both types and of diagonal links 
occupied in these configurations. 



2. Partition function and correlation functions 



With the transfer matrix T constructed, we have ac- 
cess to statistical properties of the model on a L x L T 



torus by means of the partition function, which reads 
Z {u, v,z D ) = Tr(T Lx ) = T Cl , C2 ...T Cl _ Cl 

Internal energy and specific heat - From the TM 
elements, one also gets the internal energy of the system: 

ln(T Cl , C2 ...T c ^, Cl ). 

We computed the specific heat per site c v = 
r \ d(E}/dT by a numerical differentiation of (E) with 
respect to the inverse of the coupling strength. For the 
calculations of Sec. IIV1 (E) was computed exactly for 
L = 4, 6, 8- wide systems (with L T 3> L so that c v does 
not depend on L T ). For wider systems L = 10, 12, (E) 
was estimated by searching the leading TM eigenvector 
and evaluating from it the energy per row (E)/L T (the 
consistency of this method, valid in the limit of infinitely 
long systems, was checked for L = 8 first). 

Correlation functions - Correlation functions are 
computed by means of TM iterations. Dimer-dimer cor- 
relations, for two dimers with a specified orientation and 
relative position r, arc obtained as the ratio over the par- 
tition function Z of a modified partition function Z' tak- 
ing into account only the 2D-configurations where these 
two dimers are present. Mathematically, a projector P 
onto the corresponding ID configurations is inserted in 
the matrix product in Z' at positions corresponding to 
both dimers: 

Z' _ Tr(PT r - PT L -" r -) 
W ~ ~Z ~ Tr(T^) 

where r T is the projection onto the r axis of the relative 
position r. Even if we do not discuss them in the body 
of the text, monomer- monomer correlations are also ac- 
cessible via TM calculations. In the TM language, the 
insertion of a monomer in a given row requires to replace 
the TM by a modified matrix T' which forbids the oc- 
cupation of links emerging from this site (this involves 
a change of topological sector). For both dimer-dimer 
and monomer- monomer correlations, we checked that in 
the non-interacting case, the exact results of Ref. [Ill are 
recovered. 



3. Analysis of real-space correlations 



When the system is in a critical phase, conformal trans- 
formation techniques [24| can be used to analyze the 
real-space decay of correlation functions on a L-wide in- 
finite cylinder. Consider the correlation of two dimers 



11 



(or monomers) located at a relative distance r r along 
the cylinder axis. If the exponent of the infinite-plane 
power-law correlations is 2a, correlations on the cylinder 
scale as: 



2irr T . 



C(r T ) ~ cosh( — - 1 



which can be approximated, for distances much larger 
than L/(2tt), by an exponential decay with a correlation 
length Lj (2na). We checked our calculations on the non- 
interacting case on the square lattice for a L = 6 cylinder. 
We found a correlation length £ = 1.02(2), close to the 
expected (thermodynamical limit) value L/(2wdi t o) = 
3/tt ~ 0.95. The small discrepancy can be attributed 
to the 1/r 4 correction present in the correlation func- 
tion (l4j |. Repeating the same calculations for monomer 
correlations, we extracted from the exponential decay a 
critical exponent c?o,i ~ 0.26, in good agreement (for such 
a small L value) with the exact value 1/4. 



4. CFT analysis of largest TM eigenvalues 



The CFT analysis is also useful to compute the central 
charge c, as well as the scaling dimension of the dimer di i0 
and monomer c?o,i operators, from the size dependence of 
the leading TM eigenvalues [l2j]. The central charge is 
estimated by the L-dependence of the largest eigenvalue 
Ao of T, which is directly connected to the free energy 
per site /o, in the limit of an infinitely long system: 

/o = -iMA ) = r-^ + <4). 

In practice, we add a A term to the previous fitting ex- 
pression [131 ]. which is justified by symmetry reasons and 
improves the agreement between raw data and the fitted 
expression of Jq. Estimating c thus requires to perform 



a fit of /o with at least 3 values of L. The size of the 
TM is a limitation to the number of sizes accessible; nev- 
ertheless, one can use the already mentioned topological 
and translation invariance to adress larger sizes. Indeed, 
one can assert by symmetry reasons and check numeri- 
cally that the leading eigenvalue of the TM is found in 
the q = even sector. The scaling dimensions of dimer 
and monomer exponents are determined from the largest 
eigenvalues in other symmetry sectors fl2| : 



4- 



• A5(g = 7T) - 

■Ag(g = 0). 
Ag(g = 0) N 



^ + o(^) 

2ndo i 



V- 



This allows to determine c and d e ^ n (with (e + m — 1)) 
from the power method [25| , requiring less memory than 
for a full diagonalization. To minimize finite-size effects 
on c and d e . m , wc use different estimates for these quanti- 
ties. Each estimate is the result of a fit with either 3 (for 
c) or 2 (for d e ^ m ) consecutive values of L/2. In this way 
one obtains estimates c(L — 4, L) and d e>m {L — 2, L) (for 
L < 14). The finite-size effects on these estimates be- 
come smaller when increasing L; in order to obtain more 
reliable values we make extrapolations of these estimates, 
assuming that they scale as 



c(L-4,L) ~ c* + 
d e ,m(L — 2,L) ~ d* m 



K 



K' 



{L-2f 



Such estimates and their extrapolations are displayed 
in Fig. [21 |U and \§\ In general, the c* extrapolations 
were done with estimates c(8, 12) and c(10, 14), while 
for d* e m extrapolations we used the estimates d e , m (8, 10), 
d e>m (10, 12) and <i ejm (12, 14). The difference between the 
extrapolation and the closest estimate (the one with the 
largest L) gives the order of magnitude of the error. 
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